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ABSTRACT 


A set of highly efficient computer programs based on the 
Marcum and Swerling's analysis on radar detection has been written 
in MATLAB to evaluate the probability of detection. The programs 
are based on accurate methods unlike the detectability method which 
1s based on approximation. This thesis also outlines radar 
detection theory and target models as a background. 

The goal of this effort is to provide a set of efficient 


computer programs for student usage and teacher's aid. Programs 


are designed to be user friendly and run on personal computers. 
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I. INTRODUCTION 


A. BACKGROUND 


The Marcum-Swerling (M-S) models are the most commonly 
used target models in modern radar analysis. The other target 
models have also developed over period of time, such as chi- 
Square (Weinstock), log-normal. However in this thesis only M- 
S models will be discussed. 

Histrocially, the earliest descriptions of a target were 
in terms of a single cross-sectional area value. This quantity 
was usually some type of an average cross section over the 
aspect angles which the system designer considered most 
probable. The approach led the system designers to associate 
a single cross-sectional area of unvarying value, with a 
target and thus led to the model of a steady-state target. 
Following this concept, systems were designed to achieve a 
somewhat arbitrarily specified signal-to-noise ratio at some 
specified maximum operating range. The inherent basic 
assumption was that the required detection capability and 
parameter estimation accuracy could be achieved if this goal 
was met. The dominant rule of the human observer in early 
detection systems tended to make this approach acceptable. As 


the requirements on radar systems became more demanding, and 


aS automatic processing of radar returns was developed, the 
Situation underwent a change. The desire to optimize radar 
designs established the need for more precise target models. 

The initial work of J.I.Marcum of the RAND Corporation in 
the late 1940s applied the work of Rice of the Bell Telephone 
Laboratories, relating to steady-state signals immersed in 
noise, to the problem of radar Signal detection. Marcum 
effectively gave a complete treatment to the statistical 
problem of a group of constant amplitude signal pulses in the 
presence of noise. His work resulted in an evaluation of 
previously untabulated functions and a direct application of 
the results to a gamut of problems in automatic detection. 
The statistical approach used in Marcum's studies is the basis 
for much of the work that followed in the radar detection 
area. 

The detection of signals resulting from a fluctuating 
target 1S basically different than for signals resulting from 
a non-fluctuating target. A requirement remained, therefore, 
to produce the same type of analysis and supporting tabulation 
of basic functions for a fluctuating target as had been done 
by Marcum for the nonfluctuating target. This work was done by 
Swerling, in the early 1950s. He extended Marcum's approach by 
employing four target models and two different density 
functions in conjunction with two extremes of correlation. 

Swerling's case 1 and case 2 are the target models which 


describe large complex targets, such as aircraft, rain 


clutter, and terrain clutter. They represent two extremes of 
correlation, and the statistical model is derivable from the 
Boyvstecal scattering characteristics of such bodies. It can 
also be stated that empirically derived data on such targets 
are in very good agreement with those obtained from the 
mathematical model. 

At the time Swerling was doing his work, it was realized 
that the model suitable for large complex targets would not 
give an adequate description of large simple structures. The 
problem of selecting a model to describe this type or class of 
target is a difficult one to handle directly. There is no 
Obvious parallel development for it from the scattering 
Sm@easacteristics of a specific body type. The approach 
employed was to establish qualitatively the basic statistical 
behavior of the target cross sections of interest. Having done 
this, a convenient form of distribution was selected. 

Swerling picked a class of distributions know as the chi- 
square class, which has the exponential density function 
distribution as one member. By an appropriate selection of a 
class parameter, namely, the number of degrees of freedom, the 
qualitative properties desired for the distribution associated 
with large simple structures was obtained. Swerling's cases 3 
and 4 represent targets that behave as if they have four 
degrees of freedom and are valid for targets such as rockets, 
missiles, and space-based satellites. 


Swerling's two classes of density functions evaluated at 


two extremes of correlation, together with Marcum's constant 
target model (case 5), have been the bases of virtually all 
radar detection analyses. 

Five target models according to Marcum-Swerling scheme are 
as follows: 

a) Swerling case 1 - The echo pulses received from a 
target on any one scan are of constant amplitude throughout 
the entire scan but are independent ( uncorrelated ) from scan 
to scan. This assumption ignores the effect of the antenna 
beam shape on the echo amplitude. The probability density 
function of the radar cross section oO is given by the 


exponential density function. 


w(o,G) =+el-9/8) , a>0 (1) 


Q 


where 6 is the average radar cross section 

b) Swerling case 2 - With the same density function as 
case 1 but the fluctulations are more rapid than in case 1 and 
are taken to be independent from pulse to pulse rather than 
from scan to scan. 

c) Swerling case 3 - The fluctuation is assumed to be 
independent from scan to scan, as in case 1, but the 
probability density function is given by the chi-square 
distribution with four degrees of freedom. 


—Z2 


w(o,0) = oe us (ROG = @-(Ke/e) = 49 e@-20/o (K=2} (2) 
- € 0 


d) Swerling case 4 - With the same density function but 
the fluctuation is pulse-to-pulse according to Eq(2). 

e) Marcum's model (case 5) - With the constant density 
function which represents steady-state target. 

The probability density function of equation (1), applies 
to a complex target consisting of many independent scatterers 
of approximately equal echo areas. The probability density 
function assumed in case 3 and 4 is more indicative of targets 
that can be represented as one large reflector together with 
other small reflectors. 

There are curves available that can be used to calculate 
the probability of detection for each of Swerling cases, but 
for parameters between those charted, the designer has to 
interpolate. This can be inaccurate sometimes. Therefore it 1s 
convenient and useful to provide accurate programs to 


calculate for any parameters needed by the user. 


B. RADAR RECEIVER MODEL 

dee DESCRIPTION 

The typical super heterodyne radar receiver and the 
mathematically receiver model is depicted in Fig 1. The 
difference between square-law and linear envelope detector is 
that a square-law envelope detector is used in the small 
signal optimum receiver and a linear envelope detector is used 
in the large-signal receiver. For mathematical convenience, 


the square-law detector is applied, but the performance of 


both detectors is practically the same for all signal-to-noise 


ratios. 


TYPICAL SUPER HETERODYNE RADAR RECEIVER 





Figure 1: Radar Receiver Model 

The multiple-pulse detection is used to improve the 
detectability of the signal by achieving integration gain. The 
noncoherent integration provides an integration gain even when 
the signal has a random phase and is rapidly fluctuating, such 
as Swerling case 2 and 4. By contrast, coherent integration 
needs a nonfluctuating or slowly fluctuating signal with 


predictable phase characteristics. 


2. ASSUMPTIONS 

The assumptions embodied in the M-S detection problems are 
as follows: 

e A received pulse train consisting of N samples of noise 
Or Signal-plus-noise is available. 

e The signal is imbedded in white Gaussian noise of known 
spectral density of N,/2. 

e The signal is of unknown phase type, where the RF phase 
between pulses in the train are randomly distributed. 

e The processing consists of a matched filter and square- 
law envelope detector which operates on each pulse of the 
train and a linear integration which combines the n square-law 
envelope detected pulse. 

e The directivity pattern of the antenna is rectangular so 
that the n return pulses resulting during the antenna's dwell 
time on the target are unaffected by the antenna's radiation 
pattern. 

e The target's cross section can be described by a chi- 


squared distribution with 2k degree of freedom. 
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C. RADAR DETECTION PHILOSOPHY 
Radar detection is complicated by the fact that the target 


cross section O is a random variable fluctuating with time and 


both noise and clutter. Extraction, or parameter estimation, 
1s likewise complicated by random fluctuations of the target 
echo. For the radar detection case Neyman-Pearson criterion 1S 
used which needs neither a priori probabilities nor cost 
estimates. In radar terminology whose objective is to maximize 
the probability of detection for a given probability of false 
alarm. This objective can be accomplished by uSing a 
likelihood ratio test. Specially, there exists some 
nonnegative number 7 such that hypothesis H, (1.e.,target is 
present) 1s chosen when 
Aly) = Fanl¥) , n (4) 
Lay) 

and hypothesis H, (no target present) is chosen otherwise. 
There are two types of errors which can be made. If we Say a 
target 1S present when in fact it is not, an error of the 
first kind is made (see Fig 2). That is, we choose H, given 
that H, 1S true which is denoted as the probability of false 
alarm. Similarly an error of the second kind is made when H, 


1s chosen and in fact H, is true which is probability of miss. 


Actual 


Correct decision |Missing target 


No target present |Error of second 
kind 


False alarm Correct decision 


Error of first kind |Target present 





Figure 2: Error of detection 


To optimally detect the signal with unknown phase, 


consider the narrowband signal of duration T given by 


s(t) =Aa(t)cos[w,t+6(t) +d] ce 


=S;coso-s,sing 


where S,=Aa(t)cos[wt+08(t)], S,=Aa(t)sin{wt+8(t)], Ais signal 
amplitude, a(t) is an envelope function of duration T, 8(t) is 
the signal phase modulation. Waveform r(t) is assumed to be 
prefiltered by an ideal low-pass filter that is distortionless 
within (-w.,w,) and zero outside the interval; therefore the 
filter will not affect the band-limited input signal s(t) but 


results in a band-limited statistical process. Under this 


assumption the input waveform r(t) can be described by the 


sampling theorem as 


bar SI Sei le Joly ie) i5=41,2; 3) 03 2 eee (6) 


c 
Where 2f.T 1s the number of samples, and t=1/2f,. 


The hypothesis testing can be described as: 


an 
KN 
M 


gst+n ( target present ) on 


an 
KN 
T 


n ( target absent ) 


Then under each hypothesis the probability density functions 


can be written as: 


m 3 : 2 
; » (r,-Sy) ; » Ty 
2 
aa (r) al e 20, : f (r) al e 20, 


- (/2no,)” 2 Vorname 


The likelihood ratio A(x) is given by 


exol-5* [r3-S;]*) 
Ear) : = 20 N, 
Aly) = 22 aaa 
a ae ar 
m/f£.=T fn . 
ean, 
z (9) 
exp{-+ f “[r(t)-s(t)]? de 
2 No Jo 


= eh 2 
exp{ WF, ip r(t)2dt} 


exp(-= [" s*(t) dt + < fir (t) s(t) dt} 
0 0 


where E = [eee is the energy of the signal. Substituting 


Eq(5) ainto Eq(9) yield the following expression for the 


10 


likelihood ratio A(x | ) 


A (rl) = ajo f GLE) IS 6) al) Ccos® 
No 0 








No 
(10) 
m2 , 
. (f° x(t) s(t) de] sing |} 
Let the signal-to-noise ratio % =E/N, and 
2 z 
Vat) = UGS nCG) ac 
i ae - 
(19) 
2 a 
ve) = fe CejSss) ae 
O ra Q 
uSing Eq(1l), Eq(10) can be written as follows 
A(rld) = e¥/? exp {VPly,(T) cosh - yp(T) sind]} 
(12) 


e-¥/2 exp {VP r(T)cos(p+a)} 


where O=tan'(y)/y,),7(T)=[yo(T) +y;(T)]’”7 is the envelope of the 
radar signal out of a filter matched to the waveform (2/N, 
w’)s(t), sampled at time t. Assume 6 is a uniformly 
distributed density function U(0 2m). Averaging with respect 


to @ yields the averaging likelihood ratio as follows 


A(r) ne t/2 2 [Sevier mess (oen) ag 
27% Jo 
(13) 


=e W/2T [Vpr(T)] 


where I,(w'r(T)) is the modified Bessel function of the first 


Eek 


kind and order zero. Therefore the threshold test becomes 


Hi, 
Iy(V¥ xr(T))] 2 ev? A(x) (14) 
Hy 
where eY’’*A(r) is the operating threshold. The corresponding 


receiver implementation is shown in Fig 3 
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Figure 3: Optimum detector for radar signal of unknown phase 


D. OBJECTIVE 

The Marcum-Swerling target model is the most commonly 
used model to calculate radar range performance. The objective 
of this thesis is to review the underlying theory of radar 
detection for this model and then develop a MATLAB programs to 


compute probability of detection and maximum detection range. 
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Ii. MODELS AND MATHEMATICAL METHODS 


The procedure for the calculation of average probability 
of detection P, for Marcum and Swerling cases is as follows: 
1)Find the single pulse characteristic function, C,(p), by 
transforming the single pulse probability density function 
Esun(y) . 
feria the N pulse characteristic function, C,(p)=[C,(p) ]*. 
3)Average the N pulse characteristic function over the target 
Sietribution to find the average of C,(p)( Cy,(p)). 
4)Transform C,(p) to find the average ensemble probability 
density function f,,,fy) and f, (y) for the N pulse return. 
5)Find probability of false alarm P,;, by integrating f,(y) 
maomey, tO infinity. 
6)If P;, is given find y, by means of the mathematical 
recurSive method. 
7)Find the average probability of detection P,, by integrating 
the density function of signal plus noise f,,,(y) from the 


threshold y, to 


A. MARCUM'S (NON-FLUCTUATING) STEADY-STATE TARGET MODEL 
The non-fluctuating target model is applied to spherical 
or nearly spherical objects, such as balloons, many 


wavelengths in diameter. The target model, therefore, can be 


13 


represented by a constant-valued radar cross section. 
The Rayleigh and Rician probability density function of 
the envelope detected output for N pulses are given by 


-r? 


Loa Poet ern ae sere 


(15) 
-(r;+9) 


Ber |e ean ecg ren 
where W=E/N, is the single pulse peak signal-to-noise ratio 
for a steady target. Since the noise received in the ith 
observation interval is assumed to be statistically 
independent of the noise in every other observation interval, 
and the signals are noncoherent between the different 
observation interval, the initial phase 0, in Eq(5) in the ith 
interval is also statistically independent of 9; for all 1i#j. 


The likelihood ratio test in Eq (9) can be written as: 


Ar) 222 Eat tees Eyl Sar 82 = = Sy) -f f, (24183) 
Eo (x JO) ten Fy CF ]O) 
N 
= e-*/27, (r,/P) (16) 


i1 


N 
= eo? It I, (eau) 


From Eq(16), the test statistics is 


= 


hewmen), Stee Nea) (17) 


bh 
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For small signal, the modified Bessel function I,(x) can be 


approximated for x<l1 as: 


eae. 
ue = 1+—+—_+..... 
a (*) a 64 
(18) 
a c 
Silene bei eee eee x4 
o (X) ( rn F (x*) 
Therefore Eq(17) can be rearranged and modified as 
N 
y> rian! (19) 
T=1 


where NY’ is the operating threshold which is determined by 
specifying the false alarm probability. 

Marcum utilizes a square law detector law which allows the 
Signal-plus-noise probability density function to be expressed 
directly in the signal-to-noise ratio yw. To simplify the 


Sabculation, let 


2 
z 


N r N 20 
y=) 5 =e a: >) 


f=1 eT 





which Y is compared to a suitably modified threshold ¥Y, 
Bemetind the probability density of Y, first the 
Peewability density of f,,,(q,;) is found by Jacobain 


Beamstormation f.,,(r,) from Eq(i5) into f,,,(q,;) as follows: 





a qi*¥ 
F (r,) 


ri ai 
een (3) =F ac ar I, (f2q;) , gJ;20 (21) 


i 


Since random variable Y is the sum of N statistically 
independent random variables q,, by applying the relationship 
between the sum random variable and the individual random 
variable of the characteristic function, and the Campbell and 


Foster tables of Fourier transforms ; the C,(p) 1s 


N N po ‘ 
Cy (p) elk Ca, (p) atl e7 Shp ae d),) dq; 


i=l i=1Jd -# 
Nf? jpg, .-(ay+¥/2) 
See cca I,(/2a;8) dq; (22) 
fe 
2 ENE ee 
eo 


From Campbell and Foster Tables of Fourier transforms pair 
650.0, the probability density function of f,,,(Y) can Bbegiemme 


s#n 


by taking the inverse transform of C,(p). 


fon(¥) = (22) ew? Ty, (YZNPY), 20 (23) 


For small Y, the modified bessel function I,.,(Y) can be 
approached as: 


ne (nis ae fer (24) 
2™m! 22(m+1) 2°(m+1) (m+2) 


Therefore the resulting normalized square-law detected 
probability density function is 


yN-le -Y-Ng/2 


fan ¥) = Tay 


, Yeo (25) 
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for noise only the probability of density function f,(y) is 


7 YN-lo-f 
Eee FEN =a Y20 (26) 


A computer simulation of f, (y) at the square law detector 


output is given in Figure 4: 


& 
9 
= 
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oO 
Cee 
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~~ 
= 
3 
© 
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Figure 4: Simulation of noise at suare law detector output 


After square-law detection, the normalized square-law 
detected variate (y) is summed over the N-pulses in the pulse 
train and then compared against a normalized threshold voltage 
(y,) to determine the presence or absence of a target return. 


The probability of false alarm can be obtained by 


oy 


integrating frem the threshold) (y,) Jcosimeie, 
ao N=) =v 
Pe a ae ay (27) 
yb (N-1)! 
The incomplete Pearson gamma function which is very useful 


for describing M-S model is defined as: 
+ “YyP 
T(u,p) =fYP') St ay (28) 
0 Pp! 


And in term of the Pearson's incomplete gamma function, Eq(27) 
is 
Pe UN, Vagina diay Neer (29) 


where u=(y,/N'”) and P = N-1. 
Eq(29) can be successively integrated by parts and given 


ao. 


= Yb SN ea ee 


Pra (Ne ty) = “tumayy ye Ys 


(30) 


N-1 K e*b 


Yb 
> K! 


The above equation also can be approximated for N>>1l1 by 


letting 


N-1 ,. (N-1) (N-2) i 
1S eee (= N >1 
ee ae ls ,-Nz1 (31) 
Yp 
and by applying Stirling's approximation 
N 
N! -/20N( =) (32) 


The probability of false alarm in Eq(30) can be approximated 
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as 


\ Ae ees 
pemetiee ( \yo ee NST (33) 
fa V On ( ny? Y-b-N+l is 


The probability of detection for steady-target (case 5) 
Eameoe Obtain by integrating f,,,({Y) from y, to °° and is given 


by 


S+n 


P,(Y) em f,.,(¥) dY 


(34) 
- 1- fee) (WI)/2 e-Y(W)/27,  (VENPY) 20 
The Q function is defined as 
eo N= 2 Z 
Oy (a, B) =f, a) exp[-(4—"—)] Ty. (av) dv (35) 


For N pulses the probability of detection P, can be 


Wiemeren in terms of QO function as follows: 


P,(y) =O(VZNG, /Z¥,) (36) 
Another approximation can be made by using Gram-Charlier 
series expansion ( Appendix A ) and noting that Gaussian 


family is closed under linear operation. 


eG) 
my, = (-F) 3 |p-o = N(1+p/2). 
2 
m, = (-j)2 OP) ea ey/2)2+n (ry) ae 
ape 


o?=m,-m; = N(1+W) 


where m; is the ith moment of the distribution of Y and o is 


ig 


the variance of Y. 


eae 1 @ 7 {¥-N(1+y/2)]? 
s+n ORN(L+) ee 2N(1t+w) (38) 
CY ie eee tea 2 

J2nN 


Therefore the probability of false alarm can be represented 











as: 
7 ATEN SE - =e Y= 
Pepa te WW ay = [ —e 2dzZ where Z-=—= 
Ys J2TN ,-N vom JN 
(39) 
JN 
Yous 
= o[—-—] 
JN 


Y, can be solved from the above equation 


Y, = JN O° (P,,) +N (40) 
The probability of detection can be obtained as 


(ye Setar) 


] (41) 
VN(i+¥P) 


Pp = ol 


B. SWERLING'S (FLUCTUATING) TARGET MODELS 

1. SWERLING CASES 1 AND 2 

The radar cross section of a target (0) is the area of a 
hypothetical reflector that scatters all radar beam energy it 
intercepts omnidirectionally, such that it produces an echo at 
the radar antenna equal to that energy actually received from 


the target; that is 
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J (power received by antenna) 


(incident power density per4n raa) 
(42) 


z 
slim 4nR? eel ER 
ae ED 


Where £,1s the electric field incident on the target, E, is 
the reflected field as meaSured at the radar receiver antenna, 
and Ris the distance from the target to the radar receiver 
antenna. 

The total electric field received from a complex target 


can be expressed as a summation, namely, 


te j4and 
E,=| y ER, exp (—+—) | (43) 
K=1 


and therefore the radar cross section is 


Fhe 


0-|5. TH exp (———*) ? 


K=n 


“1D (Td [cos 


(44) 








_ kK) ect BK 2 


where 


G., is the crossing section of individual scattering 
elements, 

n the number of scattering elements, 

d, the range of the Kth element,and 

X the wavelength. 

In order to relate the proceeding to a real target two 


assumptions shall be made: 


Za: 


(1) On a short term basis 4nd,/A is a random 
variable which can take on any value from 0 to 2m with equal 
Probaba live, 

(2) The individual scattering elements have equal 
radar cross sections, that 1s, Ok=6,. 

Therefore the random variables can be expressed as the x 
component (0,)'"'/*)cos (4nd,/A) and the y component (0,)7” 
cos(4nd,/A), the problem of determining the probability 
density function of oO is identical to the problem of 
determining the distance moved from an origin in a two 
dimension walk problem of n steps, where the length of each 
step is (6,)?” and the direction of each step is perfectly 
eahielomy, 

The probability of having a component (x,y) after n 
steps where n is a comparativity large number is given by 


exp [—(x-*+y-)//n0,] 


(45) 
TNO, EBOY, 


w(x,y) ax dy = 


Converting from rectangular coordinates to polar coordinates 


yields 





w(R,@) dR dO = ~ exp[-(R?)/no,] dR a0 (46) 


0 
The marginal distribution of R, obtained by integrating the 


preceeding with respect to 08,is 


w(p) aR = = Gx (ay Glee (47) 


0 
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Now, Since O=R’=x°+y* and do=2R dR Eq. (47) becomes 


_ exp[-a/na,] 
nd, 


w(a) for a> 0 (48) 


w(O0)=0 for os0. This demonstrates that the probability 
density function of the cross section is exponentially 
distributed with an average cross section Sd=no,. Thus, the 
average total echoing area is the sum of the individual 


echoing areas of the individual elements. 


a. Case 1 Scan-to-scan fluctuations of exponential 
density function. 

The scan-to-scan fluctuations of an exponential nature can 
be applied to targets such as jet aircraft when a radar having 
a fairly high pulse repetition rate and scan rate is employed. 

For case 1, Swerling derived the probability of detection 
by initially defining a target model where the received signal 


power 1s exponentially distributed, namely, 


2 exD = 7/0) 50 
w(X,X) a aaa for x (49) 


= 0 for x <0 


where x = input signal-to-noise ratio. 


Ul 


X = average of x over all targets fluctuations. 
The characteristic equation for the probability density 
function resulting from the integration of N pulse returns 


from an exponential fluctuating target, when there is complete 


es 


correlation between pulses, is derived from the non- 
fluctuating case (Swerling case 5) characteristic function as 


follows: 


NX.) 


-Nx P+1 





el ae 

-Nx(—2. ) 
es é +1 (50) 
Cy (P) ai w(x, X) oc ee 


1 
(pel) eee Ie phsa al 
For noise only, the characteristic equation is the same as 


for case 5. 


1 
(a 51) 
Y (1+p) * 


The probability density function f,(y) obtained by using 
Campbell and Foster 431, p.44, is the same as the probability 


density function of the nonfluctuating case. i1,e, 


= ese 52 
fn Ye eran (52) 


For signal plus noise, the probability density function 
f,,,(Y¥) was obtained by using Campbell and Foster pair 581.7, 


p.64. For N>l, the density function f,,,(Y) is 
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ee ely G 
o - (p+1)** [1+p(1+Nx) ] = 





Se eee ee [T (w-1) -P'(w-1, —NX_y) 
Nx P(w-1) (NK) ** NX 
L+NXx 
Pees 53 
7. 1 owe, P(N-1, —=y) eu 
ee (lt) (8 = Ne) 
eee Teeny) | (ee ND) 
Nx NX y+ JN-1 
NX 


where I'(v,z) is the incomplete gamma function and I(u,p)is the 
incomplete Pearson gamma function. 
For N=1, the probability density function can be obtained 


from Campbell and Foster pair 438, p.45 as 


meee ext) (54) 


g+n 
therefore the probability of detection can be obtained by 
integrating Eq(54) from 0 to Y, and the result is given by 


Swerling as: 


y 
1-Pp=f f..,(y) dv=I{——2— ,N-2)-(1+nX) £(¥,) (55) 


— 


From above, P,can be obtained as 





Y,, 4... N-1 
= = f N-2 1 ee 
pet en 6 )] + (1+— =) 
(56) 
(1+ (1/Nx) J/N-1] 


Zo 


Another approximation can be made by inspecting Eq(56) and 
let Nx>>1l and P,,<<1l such that the Pearson incomplete gamma 


functions are close to unity. Then Eq(56) can be rewritten 


approximately as 





al ' N-1 Nx>1 


Peat 
D ( Nx/2 Pra<i 


(57) 


Taking the logarithm of Eq(57) and using a series expansion of 


each term, 


ine N= em (bene - 





NX (NX) (1+) 
NX 
e(w-1)[—2--2__2_ 32) 58) 
(Nx) 2 (Nx/2)? 3 (Nx/2)? 
pk eee ee 
Nx/2  (NxX/2) (Nx/2) 
Taking only the first terms and using Eq(40) yields 
l ZZ 
nN Dp=-——(Y,-N+1) 
NX 
(59) 


2-4 [YN 67 (P,,) +1] 
NX 


b. Case 2 Pulse-to-pulse fluctuations of exponential 

density function. 
Pulse-to-pulse fluctuations of an exponential nature apply 
to targets such as propeller-driven aircraft (1£ the 
propellers contribute a significant portions of the echo 


area), to targets where small changes in orientation would 
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establish significant changes in echoing area (such as long 
thin subjected to a high-frequency signal), or to targets 
viewed by radars with sufficiently low repetition rates. 

For case two, the received signal return still belongs to 
Sespemential distribution. For this case the signal is 
completely decorrelated (pulse to pulse fluctuations exist), 
the probability density function for noise only (signal-to- 
noise ratio equals to 0) still the same with Case 1. The 
probability of density function of signal plus noise for case 
2 can be obtained with the same procedures, but have the pulse 
to pulse integration. The characteristic equation for the 
Single pulse is obtained by letting the characteristic 


equation of case 1 to be N=1 


= 1 
C. = ———— (60) 
1D) [1+p(1+x) ] 


For N pulses which is completely decorrelated (independent), 
the characteristic equation is just the Nth power of single 


pulse 


Sao wie io = ie 


[| _—_ -_-—_- 61 
1+p(1+x) ce 


Therefore the corresponding probability of density 


function obtained by using Campbell and Foster "Tables of 


Fourier transforms" pair 431, p.44 yield: 
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oo (N=-2) 6 (-y/ (1 
Cn =| f——_1+ _] Ne 2nPvqy ea a va 





——————— 62) 
= Tep(1+%) (ee "aT 
The probability of detection can be expressed as: 
pee tf} fey eae 
Y> 
_ iL of Se Ee ee ery) nrgy 3 LY 
== ey ae 
(N-1) 4 1+x 1+X Co 
== 7 N-1) ] 


seen ear 
(1+x) /N 


Another approximation can be made by applying that the 
Gaussian family is closed under linear operation and by using 


Gram-Charlier series expansion (Appendix A) so that 


. aGyD) a 
m, = aay) ap |e = N(1+x) ‘ 
ol (Gay) (64) 
_ (_+)\2 Y - 2 
eee) Gow lees = N(N+1) (1+x)?. 
Variance : o? = m,-mj=N(1+X)? 
The approximate probability density function then is 
-[¥-N(1+x) ]? 
fear ( a e@ 2N(1+x ° 
f27N(1+X) (65) 
£,,(¥) = Lg (-¥-M)?/2N 


Jann 


Therefore the probability of detection and false alarm can 


be obtained by integrating from the threshold level y, to » 


and are represented as: 
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eT (¥,-N(1+X) ) (Y,-N) 


ee o (eG) 
JN(1+x) | a9 | JN | 


2. SWERLING'S CASE 3 AND 4 


For Swerling's case 3 and case 4, the radar cross section 


can be described by chi-square distributions with four degrees 
of freedom. The density function is commonly associated with 


tabilized missile tankage and can be expressed as: 


x>0 (67) 


a. Case 3 Scan-to-scan fluctuations with a chi-square 


density function with four degrees of freedom. 
The characteristic equation for case 3 which represents 


the condition of complete correlation (scan-to-scan, no pulse- 


to-pulse fluctuation) is given by 


2 (pele 
a a = ee ee. 68 
ee [1+p(1+&)]" — 
The characteristic equation and the probability of density 
function for noise only remain the same as in case 1 and 2 
a 1 
Cy(p) = ———— 
niP (pei) 


(69) 
Se (SL Gen py gS eee” 
Bay) fc ane e dp Noy 


For signal plus noise the probability of density function 


Zo 





can be expressed as: 


fey, eee pe ee eJ2™PY dp 
eee 
(70) 
pel (Oy) 
(N-1)! [(1+(Nx/2) ]? 1+ (2/NX) 


where the density function comes from Fourier transform pair 


581.1 in the Campbell and Foster tables given as 


1 See ye-l e~ (pa) ¥/2 x 


(stp)“**(sta)** I(2a) (p-a)* (71) 
Me Np On| 3. ¥S0 


where 


fe 
ane 7 Fy (Pov, 2441; 2) eee 


Eq(72) can be simplified using the relationships for the 
confluent hypergeometric function, that is 


een 2) = oben (2 IN, oy tN 2 


(73) 
Zi Nee 
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Swerling uses two identities relating to confluent 
hypergeometric functions in order to expand the proceeding 
into more familiar Pearson's form of the incomplete gamma 
BUmetion, I(u,p) which is defined in Eq (30), thus Eq (73) 


becomes 
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£.,,(y) = Li+(2/Nx) 1 Ty I y N-2] @ (-¥/1+(Nx/2)) 


[1+ (Nx/2) ]? [1+2/ (Nx) ] /N=1_ 
-{N-2) [a +(2/(NeVI 7p OY 
[as CN 2) ae [1+ (2/Nx) ] J/N-1 
e~¥/ (14+NK) 4 Ve 


[tee 2 ea) 


The probability of detection can be obtained by 
integrating the density function from the threshold Y, to ~ 


as: 


=V5 N-2 m = ~2Y¥p 
Py Yp*e"? 2¥, eae + eo t> , ( Nx+2)"? @ Nx+2 
(N-2) ! Nx+2 m=0 m! NX 





2Y, 7” (75) 
2 (N-2) Set ae 
(= See Soa a -¥ ¢ Fe _2*NK_)) ; a2 
NX NxX+2 oO m! 


b. Case 4 pulse-to-pulse fluctuations with a chi- 
square density function with four degrees of 
freedom 

With no correlation, the characteristic equation for a 


Single hit is obtained by letting N=1 in Eq(52) 


C(p) = ——+"P__. 


(a+p(a+) 1" (76) 


The characteristic equation for the sum of N pulse can be 


obtained by the Nth power of a single hit 


byl 


— — N 
ee nen = ee (77) 
Ppl tts 2). | 
The inverse Fourier transform of above equation is 
obtained from Campbell and Foster transform Pair 581.1,p.64, 
this yields the following expression for the probability of 


density function in the condition of signal plus noise 





eax (2) 1+p) % 
Leen (¥) © * — 
"p [Aol 2) ) 
- (78) 
ee eee Pr) F,(-N N: =x 
(1+) 24 ae 


The confluent hypergeometric function F,(-N,N;a) can be 


expanded as: 


Pee eee) eee ene ete Nt ia, 
2! 


oy Ns 9) S14 2) 2 ste — 
ea Nie) 1 N 1! N(N+1) N(N+1) (N+2) Be (7) 


 (-1) IN! /(N-K)!]_ ak 


fe, [(N+K-1) !/(N-1) !] 


From Eq(78) and Eq(79), the density function f,,,(Y) can be 


rewritten as 


ey a ae 
3 RE ENS x/2_,* Vie (80) 
on (a4/2)?% fet tx/20 T(N*K-1) § (NRK)! KY] 


The probability of detection is obtained by integrating Eq(80) 
from the desired threshold Y, to ». From the definition of 


incomplete gamma function, the probability of detection can be 


Su 


written as: 


Po =f faen(¥) ay 
(81) 
ee T[——_2 | N+K-1] 
oe N! yet (1+x/2) JN+K 
= oat ee 
(1+x/2)%fg 2 K! (N-K)! 


With Eq(40), the probability of detection can be approximated 


as: 


VN >> (Pe,) +N N+ 1] 


K- 
aa (5) Kk” '(@va/2) JiR’ (82) 
(>) 


P 
3 ae “(1+0/2) 9& K! (N-K)! 


C. SEARCH RADAR DETECTION RANGE CALCULATION 

The Marcum-Swerling theory represented by extensive sets 
of curves from the computer programs can be used to determine 
the detection range of a practical radar by introducing 
detection loss and others parameters. There are many types of 
detection losses which have been identified so far, and when 
these are considered, reasonable predictions of radar 
performance can be obtained. 

The radar's detection range can be determined by applying 
the desired signal-to-noise ratio determined from the Marcum- 


Swerling theory and the calculated detection loss, as given by 


a3 


P.G,.G,A*G 


- 4 
Res S |e (m) (83) 
(42) °KTBFLL(—) 
N 
where 
Ruax = The maximum detection range in meters 
+ = Peak transmitter power in Watts 


X = Wavelength in meters 
G 


tx = Transmitter and receiver antenna power gains 


Q 
ll 


Average radar cross section in square meters 


Ps 
i 


Boltzmann's constant= 1.38 x 10°% J/deg. 

T = Effective system input noise temperature in degrees 
Kelvin (° K) 

B = Receiver bandwidth in Hertz. 

L = Detection system power loss factor 

F,= The receiver noise figure 

S/N= The smallest output signal-to-noise ratio 

When n pulse are integrated previous equation can be 


written as 
P.G.G,A*on 


}23 (m) (84) 


el S 
(An) Ripe es (— 

N 

where the parameters are the same as that of Eq (83) except 
that (S/N), is the signal-to-noise ratio of one of the n 


pulses that are integrated to produce the required probability 


of detection for a specified probability of false alarm. 


34 


III. SOFTWARE DEVELOPMENT 


This chapter describes the development of MATLAB programs 
for the efficient and accurate computation of probability of 
detection based on Marcum and Swerling theory of radar 
detection. The MATLAB source code is given in Appendix B, and 
complete programs are available from Professor G.S. Gill, Code 


EC/G1l, Naval Postgraduate School, Monterey, CA 93943. 


A. PROGRAM STRUCTURE 

The overall program structure is shown in Figure (5). The 
structure is that of a main menu program which calls various 
submenu programs (mscurve.m number.m number.m ) as required. 
The submenu programs, when called, will then display the 
purpose of the subprograms. When called, the subprograms will 
call the function programs to do the actual computation. It is 
possible to exit the process from either the main program or 
from the subprograms or the function programs. The advantage 
of this format is that if the user wants to change one of the 
subprograms or function programs and wants to add other 
programs, this system can accommodate it. For each subprogram 
there 1S an error detect prevention and data entry double 


check function to inform the user and restart the process. 
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( MAIN MENU ] 
MS.M 
[ SUBMENU.M ] [ SUBMENU ( SUBMENU] 
MSCURVE.M NUMBER. RMENU.M 


[ SUBPROGRAM J { SUBPROGRAM ] { SUBPROGRAM ] 
RANGE1.M 


RADAR.M POINT.M RANGE3.M 


RADAR1.M RADTHRE.M RANGE4.M 
DETECT.M 


MEYER.M 
DETFRACT.M 


[ FUNCTION PROGRAMS ] 


SWERL1.M 
SWERL2.M 
SWERL3.M 
SWERL4.M 
SWERLS.M 


[ FUNCTION PROGRAMS ] 


THRESH.M 

THRESHM.M 
PROB.M 
MARCUM.M 





Figure 5: Program structure 


1. User's guide and instruction 


To use these programs, the following MATLAB files have 
to be copied to the user subdirectory. A brief explanation of 


each is also given. 


ems.m is the main menu program and gives brief descriptions of 


the M-S model and displays the main menu. 


*mscurve.m ,rmenu.m, and number.m are the submenu programs. 


These give the purpose of the programs when called. 
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eradar.m, radral.m, meyer.m, are the subprogram used to 
integrate swerll.m, swerl2.m, swrl3.m, swerl4.m, swerl5.m 
function programs to calculate and plot the curves p,vs S/N, 
Pa VS py, and pz vs N, respectively. 
erangel.m, range3.m, range4.m, are the subprograms responsible 
for integrating and calculating the detection range. The user 
inputs the specified detection loss from detect.m. The 
rangel.m, range3.m, range4.m, will load the data from detect.m 
and calculate automatically and display the detection curves 
on the screen. 
*epoint.m radthre.m, are the subprograms responsible for 
calculating and displaying the numerical results from the 
function program radpoint.m and marcum.m, respectively. 
eswerll, swerl2, swerl3, swerl4, swerl5, prob.m, thresh.m, 
threshm.m, bound.m, radpoint.m, noise.m, signal.m, are all 
function programs responsible for calculating the data from 
the subprogram and then return the numerical value. 

A 386 or 486 personal computer is suggested for greater 
speed. To start the program type ms and press [enter] at the 


MATLAB prompt. 


2. Printing graphical outputs 
Graphics output from all programs will be stored as meta 
files automatically. The operator can print out the desired 


graphics from MATLAB by typing !gpp <filename> [enter]. Once 


Sy) 


the user restarts the programs, the previous meta files will 


be deleted automatically. 


3. Programs options 


Once the ms.m command is given, the first screen seen by 


the user 1s as shown in Figure 6. 


STHIS PROGRAM IS DESIGNED TO ALLOW THE STUDENT TO 
SVARY THE PARAMETERS OF THE VARIOUS SWERLING CASES 
$IN ORDER TO STUDY THE EFFECTS. 


CASE #: DESCRIPTION 


of of 


Returned pulses are of a constant amplitude 
over one scan, but are uncorrelated fron: 
scan to scan. 


Returned pulses are uncorrelated from pulse 
to pulse and correlated from scan to scan. 


Returned pulses are of a constant amplitude 
over one scan, but are uncorrelated from scan 
to scan. 


Returned pulses are uncorrelated from pulse to 
pulse and correlated from scan to scan. 


The static case with constant S/N and pulse 
amplitude 


% 
% 
2 
$ 
% 
$ 
$ 
% 
t 
% 
t 
$ 
% 
$ 
% 
% 





Figure 6: Main menu descriptions 


This screen gives a brief description of the five target 
models. After pressing the "enter" key, the user will see a 


second screen as shown at Figure 7. 
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MAIN MENU 


1) THE M-S CURVES | 
2) NUMERICAL DETECTION PROBABILITY CALCULATION 


3) RANGE DETECTION CURVES 


Select a menu number: 





Figure 7: Main menu 
At this point the operator has different choices to make 
depending upon what he/she wants to do. Once the operator 
chooses one of the items, the main menu program will transfer 


to the selected submenu in Figure 8. 


PROBABILITY OF DETECTION vs. S/N 
PROBABILITY OF DETECTION vs. Pfa 
PROBABILITY OF DETECTION vs. N 
COMPARSION OF M-S CURVES 


Select a menu number: 





Figure 8: Submenu Screen 


From the submenu the user will have another set of 
choices. He can choose the item he wants to study. After he 
chooses one of the items, the following selected screen will 
be seen in Figure 9, Figure 10, Figure 11 and Figure 12. The 


user can follow the instructions on the screen to key in the 


S2 


arguments he wants to study. Then a data input screen seen at 
Figure 13 will display the data for double check. After the 
completion of above procedure, a selected graphic output will 
appear on the screen as in Figure 14. After the graphics 
display, the user can press 'enter' to get the next screen as 
shown in Figure 15. At this point, the user can either choose 
to go back to the main menu or exit to print out the graphic 
display. A selected second submenu and third submenu are shown 
in Figure 16 and Figure 17; users can follow the same 


procedure to choose the items they want to study. 


LrkkkkkkkkkkekkeeeekekekeekeekkkeekkkekekeeeekeeKeKRKKKKKKKKKKKK KKK 


% THIS PROGRAM RETURNS THE PLOTS FOR THE NUMBER 

% OF PULSES AND SWERLING CASE SPECIFIED IN THE PARAMETERS. 
% THE PLOTS WILL BE STORED IN METAFILES UNDER THE NAME 

%* "RADAR.MET" FOR AN EASY PRINT OUT. 

% 


$ (A) the swerling case number has to be determined now 
Bde te HK KKK KKK KKK KKK IAEA AAKAAK KKK KKK KR KK KK 


echo off 


Enter the case number you want to study 





Figure 9: Subprogram Descriptions Screen (a) 


kkkakkkekkkekkkkkkekekkekkekekaKeKkkekeeee KKK KKK KKK KKKK KKK 


(B). The number of radar pulses the program is to 
integrate needs to be an integer between 1 and 600 
Hk ka RK KKK KKK KKK KKK K KKK REE REA KK 


cho off 


Number of Pulses to be integrated is n 





Figure 10: Subprogram Descriptions Screen (b) 
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He te te Hee te Ke KK KKK KKK KKK KAKA KEKE KEKE EKEREKEKEKREKKEEKEKKEKKEEKEREKSES 


(Gc The probability of false alarm rate curves (pfa) to be plotted 
must now be determined. Each choice of a pfa will result ina different 
curve being plotted on the graph. You need to choose the following; 


Ls The smallest pfa curve to be plotted, pfamin = ? 
27 The largest pfa curve to be plotted, pfamax = ? 
She The step size between pfamin aand pfamax, pfastep = ? 


If you wish to plot only one curve then enter the same value for 
pfamin and pfamax. 


The suggested default step size to use is that of PFASTEP = 10, 
which is quite sufficient. It is suggested that pfamin and pfamax 


be powers of 10 as that is the normal choice. 
KRRKKKKKKKKEKKEKEKEKKEEEEKEKEKEEEEKEEEKKEKEKKKKKKKKaKKKKKKKKKEKKKKKKKKRKKKKEK: 


Figure 11: Subprogram Descriptions Screen (c) 


RHRHRKKKEHRHHRRKRKKEKREKRKKRAHKRKRKRKRKRKKRKRKRKRKRKRKRKRKKKRKKRKAKKHKKRKKRKRKKRKRKRKRHHRKkR ee 
The signal to noise ratio (S/N) in dB for which you wish to 
plot needs to be determined. The choices you need = 
to make are; 


i The smallest S/N point to be plotted, sdbmin = ? 


2: The largest S/N point to be plotted, sdbmax = ? 
Remember that S/N must be entered in GB. 


a: The stepsize between sdbmin and sdbmax = ? 
RHRRAKRHRRKRKERKRKKKHEHREKRARKRAKRARHEKKRARKRARKRARRKRKRKRAKRKHRKR thee Ae 





Figure 12: Subprogram Descriptions Screen (d) 
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% 

echo off 

The case number you is 1.00 

The number of pulses you choice are 

The max false alarm probability you choice is 1.00e-12 
The min false alarm probability you choice arel1.00e-12 


pfa stepsize is 10.00 
max signal-to-noise ratio you choice is 10.00 
min signal-to-noise ratio -10.00 

stepsize is 


PARAMETERS ARE CORRECT PRESS 1 
PARAMETERS ARE NOT CORRECT PRESS 2 





Figure 13: Input Argument Checking Screen 


MARCUM-SWERLING CURVE 


n= 10 
pfa = le-08 





Probability of detection in percent 


0 10 20 


(S/N)1, signal-to-noise ratio, DB 





Figure 14: Selected Result 
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If you want to go to the same submenu 


or 


If you want to go to the main menu 


or 


To exit this whole program 


Figure 15: Selected menu 


il 
je 


ENTER CHOICE 


lt 
NO 


ENTER CHOICE 


PRESS RETURN 


SUB MENU [ NUMERICAL CALCULATION] 


1) CALCULATION OF DETECTION PROBABILITY 
2) CALCULATION OF THRESHOLD LEVEL 





Figure 16: Selected second menu 


SUB MENU -- (THE DETECTION RANGE ANALYSIS] 


1) RANGE DETECTION CURVES WITH DIFFERENT PROBABILITY OF FALSE ALARM 
2) RANGE DETECTION CURVES COMPARSION 
3) RANGE DETECTION CURVES USE THE DETECABILITY FACTOR 





Figure 17: Selected third menu 
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B. Algorithms for M-S models 
The function programs in Figure 5 are based on the Marcum 
and Swerling detection theory. 
l. Function program algorithms 
The function programs can be executed independently as 
a normal MATLAB function. The input arguments and output data 
are listed in Table [1] at the end of this chapter. All these 
programs were implemented with the help feature of the MATLAB 
environment. Typing “help < name of the function >" will 
explain how to use them independently. 
a. Threshold computation 
The detection programs start with the calculation of 


threshold y, by applying the equation (30) in Chapter 2. 


Y N-1 yx 
EP aN iat Noe ee) (85) 
VN = «OK! 


The false-alarm probability can be represented as a 
function of Nand y,. To compute probability of detection for 
all the five cases Y, is required for a given p,, . 

Since Eq(85) is a finite power series, in order to find Y, 
for given P,,, a recursive computation method has to be used 


rewriting Eq(85) 
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N-1 k N-1 


ies -Y Yp (-Y;) 
P,.(N, Y,) = ——@ 5=P2(N-1, Y,) +e 2 
=P,,(N-1,Y,) +L (N-1) 
= Vp Y ey 
Paved, ope are = Pra(N, Yp) +e 
k=0 “*° 
=P,,(N, Y,) +L (N) oe 


Y 
=P.,(N, Y») +L(N-1) = 


From Eq (86) and (87) 


ve 
a = (88) 
To ee Ny 1) = 


For a single hit (N=1) the false alarm probability and the 


relation with the following term (N=2) are 


P,,(1,¥,)=e" , L(1)=Y,e° (89) 


Therefore this relation allow each term of the expansion 
to be based on the value of the proceeding term, therefore an 
algorithm can be formed to compute the values of the detection 
threshold y, and the number of integrated pulses N. 

The first procedure employed in the algorithm is to define 
an empirical threshold level Y, then compute P(N,y,). On the 
basis of this empirically determined value y,, an empirical 


suggested Y, is given by 
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Y,=N-J/N+2.3VL(VI+/N-1), where L=-logP,, (90) 
This value is used as the starting point to compute 
P(N, y,) then compare it against the desired value of P,;, and 
the difference between Y, and Y,,,. This value of correction ay 


can be calculated by using Newton-Raphson method by noting hat 


P ra (N, Yy) 


Pes (91) 
e-tvyN"1/ (N-1) ! 


Fae (N, Yy) 


ln 
Yer = Yyt 


The procedure is repeated until the correction magnitude 
Ay/(ayty) ais indicated that Y, 1s within a sufficient 
accuracy. A computer independent algorithm notation is given 
on Figure 18. In Figure 18, the arrow implies a specification. 
The normal execution of statements is carried out line by 
line, starting at the top, but a branch may be designated by 
an arrow which results from the execution of a statement. A 
conditional branch is denoted by a colon statement, and the 
branch is executed if the comparison condition specified on 
the arrow is satisfied. Otherwise, the next statement in the 
sequence is executed. Notice in figure 18 that the program is 
terminated when the value ay/(ay+y) is less or equal to €. 
Mie value of € can be assigned to be 10° or 107%. This 


accuracy should be sufficient for application. 
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~10g,5 (Dz,) 
n--ynt2.3/L(/L+/n-1) 
Yo 

exp (-y) 


=~ Ym*Y/M 

~- Yms+Ym 

- Yms 

= (p/Ym) *1n(P/P,,) 
- YtayY 


y 


J ee 


Figure 18: Algorithmic Program for Detection 
Thresholds, y, 





The MATLAB function for executing this process is named 
prob.m and thresh.m respectively. In Figure 18 € is the 
smallest acceptable tolerance value. ay/(ayt+y) is compared 
with €, and when it is less than or equal to €, the 
computation will stop and return the value of false alarm 


probability. In thresh.m the smallest tolerance value was set 


47 


Geube 10°, this value is sufficient for desired accuracy. 
This recursive method will be used in all Swerling cases 


to find the threshold level y, 


b. Swerling case 1 algorithm 
Following equation is used to compute probability of 


getection for case 1. 














-Y 
N-1 Y b 
P, =1-I{—2_, (N-2)] + (1+) 7( —__*_, w-2) e ow 
: VN-1 NX [1+(1/Nx) /N-1] 
=o 
N-1 Y b 
=P,,(N-1,¥,)*(1+—2) [1-P,,(N-1, —2-)]e Gen 
NX os 
NX 
(92) 


It 1s obvious to see that for N=1 the detection probability is 
e YP / iN) | The algorithm in computer independent notation which 
uses equation (92) to compute P,, is shown on Figure 19. Notice 
in Figure 18 that the detection threshold (Y,) is independent 
of the target fluctuation characteristics so that the 
algorithm given in Figure 18 is used to determine Y, for all 
target types. A MATLAB source code to compute probability of 
detection of case 1 targets is given in Appendix B. The name 


@meerenii1s file is swerll.m. 
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Enter Y,,N,X 


N * X 


— 
— 
i 
- 
— 
ae 
end 
— 
— 


P2 
-Y¥/(2 +17) 

(Lily xe 

EXE (-Yb/ (tes) 
Al * A, © (1 = P,) 
Py, + Ps 

EXP [ -yp/ (xn + 1)] 


iT) ie | 


a 





Figure 19: Algorithmic Program for Swerling 1 Target, Pp; 


c. Swerling case 2 algorithm 


Eq(63) can be modified as: 
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Py =1- N= 1) 
rae (93) 
x 
Pea (Ns Ton) 


A computer independent algorithm to implement the above 


equation is shown in Figure 20. 


va 
(ase NX) 
VPx 

exp (-X,) 


aa 
_ 
— 
| = 
a 
_ 
~_ 


> Exit Pa 





Figure 20: Algorithmic Program for Swerling 2 Target, P,, 
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d. Swerling case 3 algorithm 
From Eq(75), the equation for probability of detection for 


case 3 is as follows: 








aS (1+ 1) y) 
P.,=e ‘ er for N=1 
Gas) 
2 
Vouec Zee Nx+2 .N-2 aap (94) 
= ——— +P,, (N-1, Y,) +(——) ee x2 
(N-2)! Nx+2 NX 
a Oy, yY,Nx 
ee ey Gl Sey Nee ae )]; N22 
NX NxX+2 NX+2 


An algorithm for above equation is presented in computer 


independent notation in Figure 21. 


aa 


SE Es 
te ten 


Enter Y,.N,X 


EY Tt <i iet at fi 


Vt eh FT 


t 


T 





, tart t 1 


Figure 21: Algorithmic 


N*X 

Y, 

a 

0 

Be (i) 
YM 

0 

M+ il 


ii YM 

YMS + YM 
YMS 

A 

iv 

P2 

Y (1 +2/X,) 
(i 27x)" 3 


EXP (-2Y,/ (Xy+2)) 


De 2 (N=) 7 (Kee ZY (Ke + 2) 
Al * A2 * A3 * (1 - P2) 

1 

N 

N - 2 

LE 

L-1l 

i 

(C*L) 

Cc 

2 eth B® (X42) 

Pl + P2 + P39 --------- S Exit “Po, 
EXP -2Yb/(XN+2) 

Z XN Yb/ (XN+2)? 

G (1 + H) ---------- >Exit P,, 


program for Swerling 3 Target, P,, 
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e. Swerling case 4 algorithm 


Swerling case 4 computer program modified expression can be 


obtained by applying the power series in Eq (30) into Eq(81) 
as 
‘G 
k eras 
[1-P(N+k, 1+x/2)] 


N! 
aot (1+x/2) i k!(N-k)! 


rol x] 
































eae 
2N-1 7 “a 1+ NX m-N tee * ie 
— TD e ens é ee (95) 
oF ‘Tn m=N m! x-0 K! (N-K) . 
(25s 
es = 1+NX wi (NX)« 
1-+—— 2 D 
pe Fy ea 
m=2N m! X=0 K! (N-K) ° 
Y \m 
Y ( Nx) Nx 
2N-1 =x 1. Nx’ 1+ 5 yn WN! (le 
== )*[ e Boy cee 
1+ = d, m! > K! (N-K) ! 
(a ee 
-(—“—) 4 NX 
- +x a enn 
+(1+4% N e 2 2 
' 
m*2N m 
, (Sey 
2N-1 ‘Te! ee 
=)" e 2 2 
m=0 m! 
(a2 2 
Y 
2n-1 7! ix? 1+X oa NX 


| 





z pai 2 N! 2 K i ! 
2, “ m! o> K! (N-K) ! ‘TeNx/2? L+Nx/2 | 





(ae 
Y 
2N-1 ( ax! 1+ 
2 be daha 2 
2 1+NX/2 m=0 m! 
nN NX 


z 2 K ) N- 
(1 TNT = ABT (oe (Es) 4) 
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with this expression the detection probability computation for 
this case therefore can be simplified according to equation 
(95) and avoid the computation of infinite power series. A 


computer independent algorithm is at Figure 22. 


Enter Y,,N,% 


x 

QO 

(27 -Gxye 2) 
Yp 

AY (xX Fe 2) 
-B 


e 
YM 


le aie: bey alla 


M 


YMS 

MK 

YM ° B/M 

SUM + YM (1-MKS) 

X, ° MK ° (2N-M) /2(M-N+1) 
MKS + MK 


= 
= 
a 
= 
e 
e 
a 
= 
= 





Figure 22: Algorithmic program for Swerling 4 Target. Dy, 
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f. Swerling case 5 algorithm 
In case 5, the expression of Marcum steady target model in 
Eq (38) can also be modified by substituting the infinite 


power series for the modified Bessel function I,(X) given by 


2 (96) 
ea ee 


xX 

ITy( xXx) =(= 

wi) oe 

and interchanging the order of summation and integration, one 
obtain the probability detection for case 5 as: 


a KC (Nec) Dy ae 
Pys=e oo a CD e =e (97) 


By interchanging the order of the summation results in a 


efficient computation representation as: 


Y, oo m-N 


N- 
Pog 6 EE oO EE (1 FP o-me (NIN (98) 
D5 m! m! K! 

m=Q ‘  m=N a K=Q 7 


The right hand side of equation (98) has two terms, the latter 
term is a infinite summation of power series, therefore a 
recursive evaluation is necessary for computing P,,. 

The method to handle this infinite power series 1s to 


separate above equation into two terms. Let 
eae (937) 


where L represent the integrated pulses where L represent the 


integrated pulses. Let L 2N, then P, can be represented as 
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N=1 m L m m-N 
Eyes By ecle 
x Yee Wage _ Nx 100) 
eee eee ey oe ( 
the error term P, can be represented as 


p= et 1-S e@-nx (Nx) *) (101) 
m=L+1 mm: K=0 


if let P, =P,, then the truncated error term will be P, 


ey Vs a) Fa 
P,= ye e ae eee) 
m=L+ih =0 (102) 


Lt1-N 


sa-S et a- yee son) 
Let L be large enough such that 


L+1-N = )Ke a. (103) 


Therefore 
Pi«7 P,S€ (104) 


An upper bound can be found to limit the desired accuracy 
and avoid the unnecessary computation. The computer 


independent algorithm is in Figure 23. 
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entel Younes 


AN — Po eed 
YM -- EXP (-Y,) 
YMS © YM 
< M = i 
N : M Tae eS oe 
YM — YM <= voy 
YMS © Ngaksy aa Gul 
M — M+ 1 
SUM << YMS 
AB — Bee) 
ABS &— AB 
YM YM ° Yb/M 
YMS ©  YMS + YM 
SUM <— SUM + YM (1 - XBS) 
M — M+ 1 
K — M-WN 
XB — XB ° XN/K 
ABS XBS + XB 
<eé : Cie io) ele bs) 
P — SUM <-------------- > eXIT WITH Pps 


Figure 23: Algorithmic program for Swerling 5 Target, P,, 


5) i) 


2. Solving for round off error 

Round off error can take place during the calculation of 
case 5. That is because during the calculation of probability 
of detection, if the signal-to-noise ratio x is zero, the 
output detection probability will become the probability of 


false alarm, then the criterion in Eq (103) will become 


2 rg (Nx) Ke -Nx 
e(L)<(1 a ar) a y L 
ee ye eee (105) 
k! 


Since we can not set the minimum acceptable error to be 
zero, the computation process will not terminate. Due to this 
reason, trade-off between the accuracy and the maximum value 
of output data should be made. An algorithm to compute the 
threshold y, to deal with the round off error is in Figure 24. 
The corresponding MATLAB file name is Marcum.m which can be 
used to compute the maximum acceptable threshold level y, when 
the value of signal-to-noise ratio is zero and the number of 


pulses is large. 
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< Enter Ps,, Yo 
x : B 
M - 0 
YLX €& 0 
YLOG — LOG (Y) 
TSUM <— EAE (ay) 
M <_ M+1 
M — N-1l 
SUNK : LOG (M) 
YLX €& YLX+YLOG~SUNK 
TSUM ¢< TSUM+EXP(YLX-Y) 









> vy <- Y+LOG(XP,,*TSUM) *TSUM/EXP (YLX-Y 
— 

|Y-y¥, | E 

ce = Y, PD ND I 9 9 


Zee Gh LOG (Y) 
Y-YLX : B 
M — M+1 
SUNK < LOG (N) 
YLX €& YLX+YLOG-SUNK 
Ula c EXP (YLX~-Y) 
Yh & Y, coc c cc ccc crore > Y, 


Figure 24: Algorithmic Program for Detection Thresholds, Y, 
with round off error preventation 


og 


3. Solving for underflow and overflow 

As mentioned before, the maximum acceptable number for 
PIE AB 21s nearly e7'°?-7827428s228sosmv", the underflow and overflow 
problem arise due to the fact that the power series in M-S 
equations has the form 

on ; en x. (106) 

If the calculation for this power series 1s over the 
maximum or below the minimum number that MATLAB can represent, 
then the MATLAB will return a string named NaN or empty matrix 
({ }) and the whole computation process will be meaningless. 
If the computation needed large integrated pulse number and 
large threshold or signal-to-noise ratio, a special effort 
has to be made to retain the accuracy. One method is let the 


series be of the form 


a: k 
ey + -exp [-y+kIn(y) =e. ln (N) } =e~6 (107) 


Nel 


and let the value e? be compared to the smallest acceptable 
value e@ 7709. 782712893384040999999 in MATLAB, of course e@ 7 709 .782712893388060999999 
is much smaller than the tolerance number € assumed in the 
program. If Bis greater than 709.7821... then increase the 
number K in above equation until it is less than 709.7821... 
and then start the summation process. This method will save 
MATLAB computation of power series with a large number and 


avoid the under flow problems. An algorithm for using this 
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method to compute the detection probability is in Figure (Zao0 


and the MATLAB source code is threshm.m and iS given in 


Appendix B. 


YM Y¥/M 

(AS + YH 

YMSs 

Ro} 

Moet 

YM ‘¥/M 

YMS + YM 

SUM * YM (1 © XRS) 
XR NX/R 

XR © XRS 

(2 - YMS)(1 - xRs) 
SUM 

o 


lny 

Mel 

La(M), 

MLY ¢ LY - LAP 
¥ - MLY 
exp(°YoMLy) 


MeoN 
Rel 
XR NX/R 


XR * XRS 

R 

Rel 

0 

in (NX) 

Rey 

LAR 

(REX « LX - LAP) 
NX © REX 

EXP (-(NX - RLX)} 
BeorR 


er 





Figure 25: Algorithmic porgram for computing the probability 
of dtection,P,, with underflow preventation 
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4. Input arguments and output data 
The table below summarizes all input and output data for 


the programs. 


INPUT DATA OUTPUT DATA 
= 
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xy 





or P,, 


RANGE. 3 Dea fe oe ae CHART R vs Pd 


Table 1 : Input arguments and output data 


Where 


N=Number of pulses to be integrated 


X =average signal-to -noise ratio 
Pd= probability of detection 
Pra=probability of false alarm 

R= the desired detection range 


ns=Swerling case number 
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IV.RESULTS 

The M-S detection probability curves can be obtain by 
applying the MATLAB programs in Chapter 3. The results can let 
the user evaluate the detection probability detection 
performance easily. The detection range curves can illustrate 
relationship between the detection probability and _ the 
detection range by input the properly detection loss. 
A. PROBABILITY OF DETECTION CURVES 

As mentioned previously, the M-S model arise from the 
different fluctuations of target cross section. MATLAB 
programs can be used to plot probability of detection versus 
per pulse signal-to-noise ratio for given probability of false 
alarm for five cases. It can also be asked to determine per 
pulse signal-to-noise ratio required for given P, and P,,. 
This required signal-to-noise ratio can be used to compute 
maximum detection range from the radar range equation. 

Fig 26 and Fig 27 compare the five target models for a 
false alarm of 10°°, and the number of integrated pulse as N=10 
and N=100 respectively. When the probability is larger than 
0.33, all four cases in which the target cross section is not 
constant requires greater signal-to-noise ratio than the 
constant cross section case. This increase in signal-to-noise- 
ratio will cause a reduction in detection range. Therefore if 


the characteristic of the target cross section are not 
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properly take into account, the actual performance of the 
radar might not measure up to the performance which is 
predicted from the constant target cross section. A greater 
Signal-to-noise ratio is required when the fluctuations are 
correlated pulse to pulse (case 1 and case 3 ) than when the 
fluctuations are uncorrelated pulse to pulse(case 2 and case 
4). Fig 27 also indicates that when the number of integrated 
pulses is large, the case 2 and case 4 will approach to the 


constant target case (case 5). 


aed 
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heed 
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(S/N)1, signal-to-noise ratio, DB 





Figure 26: Probability of detection curves for five target 
models 
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Figure 27: Probability of detection curves for five target 
models 


Since both the false alarm rate P,, and detection 
probability P, are specified by the system requirement, the 
radar designer computes required signal-to-noise ratio from M- 
S curves and uses this to compute the maximum detection range. 
The greater the number of pulses integrated, the greater is 
resulting overall signal-to-noise ratio. This results in 
greater probability of detection, but it will require longer 
scan time. Figure 28 and Figure 29 give plots of P, VS SNR 
for selected P;,'s for case 1 for N equal to 100 and 500 


respectively. 
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Figure 28: Probability of detection curves for Swerling case 
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Figure 29: Probability of eerectzon curves for Swerling case 
1 





Figure 30 and Figure 31 are the Swerling case 2 plots of 
P, Vs SNR with 10 and 100 pulses integrated respectively. From 
Figure 30, signal-to-noise ratio of 4.8 db approximately is 
required to yield a probability of detection of 0.8 with 10 


pulses integrated and probability of false alarm of 10°°. 
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Figure 30: Probability of detection curves for Swerling case 
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Figure 31: Probability of detection curves for Swerling case 
2 


Figure 32 and Figure 33 show the Swerling case 3 with 10 


and 100 pulses integrated respectively. 
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Figure 32: Probability of detection curves for Swerling case 
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Figure 33: Probability of detection curves for Swerling case 
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Figure 34 and Figure 35 show Swerling case 4 pulse to 


burse, stluctuation with 10 and 100 pulses 


respectively. 
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Figure 34: Probability of detection curves for Swerling case 
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Figure 35: Probability of detection curves for Swerling case 
4 val 


Figure 36 and Figure 37 are Swerling case 5 steady-state 


target with 10 and 100 pulses integrated respectively. 
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Figure 36: Probability of detection curves for Swerling case 
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Figure 37: Probability of detection curves for Swerling case 
4 a2 


B. THE DETECTIOON RANGE CURVES. 

The ASR-9 radar is used as a sample to compute maximum 
detection range. ASR-9 radar is designed for surveillance of 
the terminal area around airports. A significant feature of 
this S-band radar is its use of moving target detection (MTD) 
processor which provides coherent integration among other 
features. The parameters of the ASR-9 radar are given in table 
below. Also tabulated is the empirically determined search 
detection range on a 1m target for the various Swerling 


models with a P, =0.9 and P,, =10°°. 


et 
ASR-9 PARAMETERS AND RANGE CALCULATION 


P.21200 kW ranges 

£=2900 MHz R, 239.41 nm 

T#1.05 ps R, #52.41 nm 

Oo #1 R,=46.68 nm 

G.=33.5 db R,=57.29 nm 

G.=33.5 dB R,=62.63 nm 

NF=#5 dB 

L#12dB 

6=1.3° Losses 

6°=75°/sec 

PRF#1200 pps Transmitter 2 dB 

Boopp™@i50 Hz Receiver 2 4B 

P20 .9 Mismatch 1 dB 

P,,=10°° Integrator 1 4B 

pal Collapsing 1 dB 
Beam Shape 3 dB 
MTI -2 AB 
Total 12 dB 


Table 2 : ASR-9 PARAMETERS AND RANGE CALCULATION 
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Figure 38: Probability of detection curves for five target 
models 
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Figure 39: Range detection curves 
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C. Numerical data output from probability of detection 


programs. 


A table of values below is included for the purpose of 


checking the programs. 


13.8155055 ey .0361810926687 . 2848035870497 -65475592794524 . 8721557722127 ee aed 


——— 
ss 
[ree | aime [mare [oman | ge 
Te Teer rrr [rms [rer | omega 


| 19.1291681 m7 .1971745188135 -5769062596383 -83647025963833 9446918645347 eae eee 









[ever [aan [aie | men | eee 
anions [eae [sani [ oman | apa 
Tiassa [sama | oranie [ann | ome 
er ere 


i 7013405 a he . 4855434894507 -7911151202538 . 9282416527967 . 9765960615308 gota eee 

















ea .73389703346785 . 9989667533193 .99999988582414 9999999999973 ae 


Em | .§6937472704103 .9139452495156 .9800436879910 .9988084719731 eS ee 
ee .78178938715194 .9999629142286 .99999999999768 ee 


{s | .8355167894356 .9999991524271 .99999915242839 .9999991524283 a’ ae 


et x13.162278 oe = 13162278 Pn 7 x7316.2278_| 


63.5481801 Poa .69852634183006 .8917070566484 .96429073448000 . 9885538115323 -99639654302453] 
a a . 8147880711125 .9759019998584 .99728588213879 .99971801636709 .99997145717056| | 
ana .999998744645916 . 9999992846098 .99999928460987 .99999928460987 9999284460987], 
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V. CONCLUSIONS AND RECOMMENDATIONS 


A. CONCLUSIONS 


Radar detection theory is outlined in this thesis and 
computer programs are developed for probability of detection 
calculations in MATLAB. These programs are based on radar 
detection theory as developed by Marcum and Swerling. These 
programs give more accurate results than the commonly 
available programs based on the detectability method which is 
empirical is native. A user's guide and instruction are 
included in the thesis. Simple examples of the programs usage 
are illustrated in the thesis. 

Students can use the programs written for this thesis to 
investigate radar system performance. These programs are cost 
effective, convenient to use and easy to reproduce since they 
are run on personal computers that are readily available to 
students. In particular, the program can aid students by 
removing a major computational burden to allow him to perform 
real world detection calculations. These programs can also be 
used for calculations in radar development work. 

B. RECOMMENDATIONS 


Due to an overflow underflow problem the programs are 
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limited to integration of 600 pulses in the probability of 
detection calculations for all 5 cases. Most of the time this 
does not present any limitation as most radars integrate fewer 
than 600 pulses. However it may be desirable to remove this 


limitation by employing a Chernoff bound or other methods. 


oy 


Appendix A 
The gram-Charlier Series is a series expansion of a 
probability density function in terms of the Gaussian 
function, its derivatives and the moments of the original 
density function. 
If P(x) is a probability density function that is nearly 
Gaussian and x has zero mean and unit variance, then the Gram- 


Charlier series expansion of P(x) is given by 


P(x)=)7 a; 6") (x) where 


120 


(A-1) 
) (x) =p (x) =e **/2 
fon (A-1.a) 
(4) _ d7o (x) 
ot") (x) ra 


Coefficient a; can be evaluated in terms of the Hermite 


polynomials: 
oi (x) =A 9-7/2) 7, (x) (A-2) 
v(2ny 
Some typical polynomial are 
Hei H, (x) =x4=6x7+3 


H, (x) =x, H. (X) =X?-10x7+15x 
oe) 5 (x) : een 


H, (x) =x? =1 7 OH, (Xt) Hxca 5 eee oes ls 


HAA X)\ Sx =x) Heya | = 2a as OI eae a 


A recursive relation for these polynomials is 
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Heese eerie) = (x) (A-4) 


A particularly useful fact is that the hermite polynomials and 


the derivatives of the Gaussian function are biorthogonal. 


Maat 1s, 
[93 (oH, (x) ax 
=(-1)? i! 4,, (A-5) 
_ fb)? at, 1=j 
0 1#¥J 
Bow evaluate coefficient a,, both side of Eq({A-l1) are 


multiplied by H;(x) and result is integrated from -» to ©. 


With Eq(A-3),the result simplifies to 


a=) f"p(x H,(x) ax (A-6) 


The first three coefficients a, a, and a, can be found by 


substituting Eq(A-3)into Eq(A-6) 


ay=[ P(x) dx=1 


a,=-f— xp (x) dx=0 (A-7) 
a= Guides dx== (1-1) =0 


«a 


With these results, the Gram-Charlier series expansion of p(x) 
in Eq(A-1) can be simplified to 


when the mean of random variable x is not zero and/or its 


i 


P(x) =O (x) +)" 7 (x) (A-8) 
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variance Q@isS not unity, the Gram-Charlier expansion of p(x) 


1s obtained by expanding a related density function g(t) ina 


Gram-Charlier series, where 


“I 


(A-9) 


I 
2p 


By a simple transformation of variables it follows that 





E 
Dix = we (ese 
dat (A-10) 
_ il -X 
ack : 
If the Gram-charlier series expansion for g(t) is 
g(t) =\° co (€) aa 
10 


then from Eq(A-10) and Eq(A-11), the expression for p(x) 


becomes 
at : i, XX (A-12) 
P(X) ae ae 


If both sides of Eq(A-12) are multiplied by H;[(x- x ) / oJ 
and the result integrated from -~ to », coefficient c, can be 


evaluated with relation (A-6) as 
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a X-X A-13 
c;- [ plo, =) dx ) 


Evaluating Eq(A-13) for the first few coefficients c; yields 








Cy=1, C= (4-3) 
Cree 0, Ge== BT} Kes LOes (A-14) 
C=O, Ci = (p71 50,450) 


(6!) 


where @ is the ith central moment of p(x), normalized by @: 


a = — fi (Sg Jol Soa (eb (A-14) 


The @ can be expressed in terms of the conventional moments 


of the distribution p(x); thus if m, denotes the nth moment 


emo x) that is 


m,=[_x*p(x) dx (A-15) 
then om through @ are related to noncentral moments m, by 


_ mM, -3m,m,+2m; 
a g2 

_ m,~4m,m, +6m,m; -3m; 
4 
"2 2 2 5 Ae) 
_ M,~IM,M, +10m,m, ~10m,m,+4m, 
a i 
ane Mm, -6mM.mM, +15m,m, ~2 Om,m; +15m,m; -5m; 

XK 

G 


In summary, to obtain a Gram-Charlier series expansion, 


the moments m, are found from Eq(A-16), or they can be 
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conveniently computed from the characteristic equation of the 
distribution. Next the ao, are evaluated using Eq(A-17). The oO, 
are then substituted into Eq(A-14) to obtain coefficient c, of 


Gram-Charlier expansion. 
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APPENDIX B 


*FILE NAME: MS.M 

Clear 

elg 

elke 

echo on 

*THIS PROGRAM IS DESIGNED TO ALLOW THE STUDENT TO 

S*VARY THE PARAMETERS OF THE VARIOUS SWERLING CASES 

%IN ORDER TO STUDY THE EFFECTS. 

% CASE #: DESCRIPTION 

% 1. Returned pulses are of a constant amplitude over one 

scan, but are uncorrelated from scan to scan. 

% 2. Returned pulses are uncorrelated from pulse to pulse 
and correlated from scan to scan. 

3. Returned pulses are of a constant amplitude over one 
scan, but are uncorrelated from scan to scan. 

4. Returned pulses are uncorrelated from pulse to pulse 
and correlated from scan to scan. 

% 5. The static case with constant S/N and pulse amplitude 
echo off 

pause 

elc 

k=menu('MAIN MENU', 'THE M-S CURVES |',.. 

‘NUMERICAL DETECTION PROBABILITY CALCULATION',.. 

‘RANGE DETECTION CURVES ' ) 

1£ k== 

mscurve 

elseif k==2 

number 

elseif k== 

rmenu 

elseif k==4 

density 

end 

% Now restart the process with some different choices 

%* Clear the workspace of unnecessary data to avoid 
conflicts 

clear 1; clear sav; clear sdb; clear vv, clear pfa, clear 

p;clear ns; 

clear 1; clear savl; clear sdbl; clear vv, clear pfal, clear 

D; 

%* Check and see if any info is to change 

elec 

echo on 

%$ If you want to restart the process ENTER CHOICE = 1 
% 

% Or 


oe =P 
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* 
%* To exit this programPRESS RETURN 
echo off 
choice = input('CHOICE = ') 
Lf eno ree -— 
ms 
else 
end 
Gic 
echo on 
% You have chosen to exit the program 
% 
% Bye ! 
echo off 


%$ FILE NAME: RADAR.M 
Cic 

format long 

echo 


Etter reeked eked a 


% THIS PROGRAM RETURNS THE FLOTS FOR THE NUMBER 

% OF PULSES AND SWERLING CASE SPECIFIED IN THE PARAMETERS. 
% THE PLOTS WILL BE STORED IN METAFILES UNDER THE NAME 

%$ “RADAR.MET" FOR AN EASY PRINT OUT. 


% 
%* (A) the swerling case number has to be determined now 


Exeter eererrerke eee eee eee ee errr kr rk rr rr er kr er ee ee 


echo off 
ns= input('Enter the case number you want to study '); 


echo on 


$5 te te tee te te tee ie ie ie ie ie ie ie te ie oie ie oie te ir ie ir ie ie ir ie ir er ier rie rr ier ie ir ie ir rie ie ie Wer er ee we oe 


sg 


toes) The number of radar pulses the program is to 
% integrate needs to be an integer between 1 and 
600. 


$F We We We te ie te ete ae ie ie ie ie ie ie ie ie ie ie ir ie te or ir ie ie ie ir ie er ee ie ro or ir ir rrr wr rire wr re ee oo 
echo off 

n=input('Number of Pulses to be integrated is n = '); 

Cle 


echo on 
Eee errr 
= eye The probability of false alarm rate curves (pfa) 


to be plotted 

%$ must now be determined. Each choice of a pfa will result 
in a different 

% curve being plotted on the graph. You need to choose the 
following; 
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ee 
ee ers = 


= ? 


% 1. The smallest pfa curve to be plotted, pfamin = 

% 

% 2. The largest pfa curve to be plotted, pfamax = ? 

% 

% 3. The step size between pfamin aand pfamax, pfastep = 
i 

% 

% If you wish to plot only one curve then enter the same 
value for 

% pfamin and pfamax. 

% 

% The suggested default step size to use is that of 
PFASTEP = 10, 

% which is quite sufficient. It is suggested that 


pfamin and pfamax | 
% be powers of 10 as that is the normal choice. 


9 ow ie te ie We ie ie ee oie te ie ie ire ie ie ie ee re ie ie ier ie or ier we re rier i wr rk er we we we we ew wo we ew 


echo off 


pfamin =input('pfamin ) 
pfamax =input( 'pfamax 
pfastep=input('pfastep = ') 


Sic 

echo on 

EK eek eke kere eee eked eee eee eee eee eee ree ee ee ee ee 
oe (D) . The signal to noise ratio (S/N) in dB for which 


you wish to . 
plot needs to be determined. The choices you need 
to.make are; 

1. The smallest S/N point to be plotted, sdbmin = ? 

2. The largest S/N point to be plotted, sdbmax = ? 


Remember that S/N must be entered in dB. 


dP dP dP DP oP dP OP CP CP OP 


a7 The stepsize between sdbmin and sdbmax = ? 


Be we wee ek ee ek i ie ee i ew ee ee ee ee ee ee oe ke ee we 


echo off 

sdbmin =input('In dB the sdbmin is ‘') 
sdbmax =input('‘In db the sdbmax is  ') 
sdbstep= input('In db the sdbstep is ' 
elg 

elec 

echo on 


) 


echo off 
fprintf('The case number you is %8.2f\n',ns); 
fprint£('The number of pulses you choice are %8.2f\n',n); 
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fprintf£('The max false alarm probability you choice is 


%8.2e\n',pfamax) ; Cio 
fprint£('The min false alarm probability you choice 


are%8.2e\n',pfamin) ; | | 
fprintf£('The pfa stepsize 15 467.2£\n yoetaseer a | 
fprintf£('The max signal-to-noise ratio you choice is 


$8.2£\n',sdbmax) ; . 
fprintf£('The min signal-to-noise ratio %8.2f£\n',sdbmin) ; 


fprintf('The S/N stepsize is %8.2f\n',sdbstep'); 
echo on 


IF THE PARAMETERS ARE CORRECT PRESS 1 
$IF THE PARAMETERS ARE NOT CORRECT PRESS 2 


echo off 
choice=input('CHOICE= '); 
if echnorece==2 

radar 

elseif choice== 

end 


echo on 
$e wwe we ee We we ie ee ee ee ie ie rie ae ee ee i oie ae ir ie rie ie ie ie ie ie wr ir ire ie rie ie ier wr ie irr we wr ir we wr we wr 


% Now run the calculations for the above data and plot it 
$ ow we we We ie We we ie ie ie ee oe ie ie ie ir ie ie ie ie ae ie ie ie ae ie ie ee ee ie ir ie ie ier ie ie ir ir ie ie ie re ir ee ie rr wr wr we wr 
echo off 
delete radar.met 
format long) 
nmax = 700; | 
axis ('square') ; | 
vv=(sdbmin sdbmax 0 100]; 
axis(vv); 
LE mol Pn 6008 Pn = tact 
error('Number of pulses input exceeds dimensions ') 
end 
if nsec ens > 5S |ense-tasctno 
error('Swerling case input does not correspond to 
allowable choices') 
end 


pfa=pfamax; 
sdb = sdbmin:sdbstep:sdbmax; 


kk=length(sdb) ; 
sav = 10 .*(sdb/10); 


We gc = 
% swerlingl 
while pfa >= pfamin 
Fores lekk 


$1 = 1:sdbmax-sdbmin+l 
p(i)=swerll(pfa,n,sav(i)); 
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= rere rien. 
a — ee 
ee 


p(1)=100*p(i); 
end 


pPuec(sdb, p, -" )-grac;tatie( ‘SWERLING 1'); 
xlabel('(S/N)1, Signal-to-noise ratio, DB'), 
Ylabel('Probability of detection in percent ') 


text (sdbmin + 1/20*sdbmax,95 ,({'n = ',int2strin)]) 
text(sdbmin + 1/20*sdbmax, 90,[('pfa = 
',mnum2str(pfamax),' to ‘',num2str(pfamin) } ) 


text(sdbmin + 1/20*sdbmax, 85,['pfastep = 
‘',mum2str(pfastep) }); 
gtext(({num2str(pfa) }); 


hold on 

pfa = pfa/pfastep; 
end 
pause 


meta radar 
elseif ns== 
$swerling2 
pfa=pfamax; 
while pfa >= pfamin 
feorei— i-kk 
% for 1 = 1:sdbmax-sdbmin+l 
p(1)=swerl2 (pfa,n,sav(i)); 
p(i)=100*p(i); 
end 
plot(sdb,p,'-'),title('SWERLING 2'); 
grid, xlabel('(S/N)1, signal-to-noise ratio, DB'), 
ylabel('Probability of detection in percent' ) 
text (sdbmin + 1/20*sdbmax,95 ,[('n = ',int2str(n)]) 
text (sdbmin + 1/20*sdbmax,90,{'pfa = ',num2str(pfamax), ' 
to.’ ,numZstr(pfamin) } ) 
text(sdbmin + 1/20*sdbmax, 85,['pfastep = 
',numZ2str(pfastep) } ) 
gtext(({num2str(pfa) J); 


ne Leeson 

pfa = pfa/pfastep; 
end 
pause; 


meta radar 
elseif ns== 
$swerling3 
fe == 
error('algorithm for swerling 3 does not work for n=1') 
end 
while pfa >= pfamin 
FOr i: kis 
for 1 = 1:sdbmax-sdbmin+l 
P(1)=swerl3 (pfa,n,sav(i)); 
p(i)=100*p(i); 
end 
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plot (sdb,p,'-'),title(' SWEREENG se 

grid, xlabel('(S/N)1, signal-to-noise ratio, DB'), 

ylabel('Probability of detection in percent' ) 

text(sdbmin + 1/20*sdbmax,95 ,[{‘n = ',int2str(n)]) 

text (sdbmin + 1/20*sdbmax,90,(['pfa = ',num2str(pfamax), ' 
to ',num2str(pfamin) }) 

text(sdbmin + 1/20*sdbmax,85,{'pfastep = 
',num2str(pfastep) ] ) 

gtext((num2str(pfa)]); 


nNoldon 

pfa = pfa/pfastep; 
end 
pause; 


meta radar 
elseif ns== 
$swerling4 
while pfa >= pfamin 
EOrs t= kk 
for i = 1:sdbmax-sdbmin+l 
p(1i)=swerl4 (pfa,n,sav(i)); 
D(i)= 2004 p17) 
end 
plot(sdb,p,'-');title('SWERLING 4'); 
grid; xlabel('(S/N)1, signal-to-noise ratio, DB'); 
Ylabel('Probability of detection in percent') 
text(sdbmin + 1/20*sdbmax,95 ,{'n = ',int2str(n)])) 
text (sdbmin + 1/20*sdbmax,90,['pfa = ',num2str(pfamax), ' 
to ',num2str(pfamin) }) 
text (sdbmin + 1/20*sdbmax, 85, {'pfastep = 
',num2str(pfastep) J) 
gtext ({num2str(pfa)]); 
hold on 
pfa = pfa/pfastep; 
end 
pause 
meta radar 
else 
$swerlingS 
while pfa >= pfamin 
for t2=1-kk 
% for 1 = 1:sdbmax-sdbmin+l 
p(1i)=swerl5 (pfa,n,sav(i)); 
D(i)=100*5 (3) 
end 
plot(sdb,p,'-'),title('SWERLING 5'); 
grid, xlabel('(S/N)1, signal-to-noise ratio, DB'), 
ylabel('Probability of detection in percent') 
text (sdbmin + 1/20*sdbmax,95 ,({'n = ',int2str(n)]) 
text(sdbmin + 1/20*sdbmax,90,('pfa = ',num2str(pfamax) , | 
to ',num2str(pfamin) }) ) 
text (sdbmin + 1/20*sdbmax, 85,['pfastep = | 
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', numZz2str(pfastep) ) ) 
gtext (([num2str(pfa)]): 
Jo(etie aja 
pfa = pfa/pfastep; 
end 
pause 
meta radar 
end 
Hota off 


SFILE NAME: SWERL1.M 

function pd=swerll(pfa,n,sbar) 

Exwrneeekeeereerrerkewkkekeee ere reer rer kr rrr ree ee ee © 
%* RETURNS THE VALUE OF THE PROBABILITY OF DETECTION FOR THE 
SWERLING1 

% CASE, GIVEN THE NUMBER OF PULSES n, THE PROBABILITY OF 
FALSE ALARM 

%* PFA AND THE AVERAGE SIGNAL TO NOISE RATIO S/N. 

% EXAMPLE: SWERL1(PFA,N,SBAR) ; 

EX KKK KKK reek eee keke keke keke kere eee ee ee ee ee 
format long 

yb=thresh(pfa,n) ; 

yY=yb/ (l+n*sbar) ; 

cte=(l+1/n/sbar) ; 
pd=prob(n-1,yb)+cte* (n-1) *exp(-y) *(l-prob(n-1,yb/cte) ); 
return ‘ 

end 


% FILE NAME: SWERL2.M 

I hk, Go ce ROR Seen aan 
% RETURNS THE VALUE OF THE PROBABILITY OF DETECTION FOR THE 
SWERLING2 

% CASE, GIVEN THE NUMBER OF PULSES n, THE PROBABILITY OF 
FALSE ALARM 

% PFA AND THE AVERAGE SIGNAL TO NOISE RATIO S/N. 

% EXAMPLE: SWERL2(PFA,N,SBAR) ; 

Err eeeeekeereeer keke eek kkk keke k kee kk rrr kr ee we 
function pd=swerl2 (pfa,n,sbar) 

yb=thresh(pfa,n) ; 

y=yb/ (1+sbar) ; 

pd=probi(n,y); 

Becurn 


% FILE NAME: SWERL3.M 

ERK KKK KK hhh kee KKK KKK KKK KKK Khe KKK Keke kkk kk 
% RETURNS THE VALUE OF THE PROBABILITY OF DETECTION FOR THE 
SWERLING3 

% CASE, GIVEN THE NUMBER OF PULSES n, THE PROBABILITY OF 
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FALSE ALARM 

% PFA AND THE AVERAGE SIGNAL TO NOISE RATIO S/N. 

% EXAMPLE: SWERL3(PFA,N,SBAR) ; 

Gere er were eee eee ee eee ee eee eee ee ee ee eee 
EuNCtiOon pe=Swenls (mia, om seam, 

yb=thresh(pfa,n) ; 

y=l1; 

for k=l:l:n-2; 

y=yb/k*y; 

end; 

yl=2*yb/ (n*sbar+2) ; 

y=y*exp(-yb) *yl; 

y2=((n*sbar+2) /n/sbar) *(n-2) *exp(-yl) *(1-2* (n-2) /n/sbar+yl) ; 
y2=y2*(l-prob(n-1,yb*n*sbar/ (n*sbar+2))); 
pd=y+prob(n-1,yb)+y2; 

recurna 


% FILE NAME: SWERL4.M 


Ber wkekekkkkeekkekkkkkekekekeKeKkkkek kkk kkkkkkkkk kerk eK eK ee © 


% RETURNS THE VALUE OF THE PROBABILITY OF DETECTION FOR THE 
SWERLING4 R 
% CASE, GIVEN THE NUMBER OF PULSES n, THE PROBABILITY OF 
FALSE ALARM 
% PFA AND THE AVERAGE SIGNAL TO NOISE RATIO S/N. 
% EXAMPLE: SWERL4(PFA,N,SBAR) ; 
EeKKK KKK KKK Kerk kkk kk kk kkk KKK KKK KKK KKK KKK KKK KKK KKK KKK KKK Kae Ke 
function pd=swerl4(pfa,n, sbar) 
yb=thresh(pfa,n) ; 
v=e2*yb/ (sbar+2) ; 
ZK=(2/ (sbar+2))%“n; 
YM=exp (-v) ; 
YMS=YM; 
for M=l:n-l 
YM=YM*v/M; 
YMS=YMS + YM; 
end 
SUM=YMS; 
OKS=ZK; 
For tien (2 m=) 
YM=YM*v/M; 
SUM = SUM + YM* (1-ZKS) ; 
ZK = sbar*ZK* (2*n-M)/(2* (M-n+l)); 
ZKS = ZKS + ZK; 
end 
pd=SUM; 
return; 


FILE NAME: SWERLS .M 


EeKKK KKK kkk kerk keke kkk keke kkk keke kkk kkk Keke eee eke 


% RETURNS THE VALUE OF THE PROBABILITY OF DETECTION FOR THE 
SWERLINGS 


| 


% CASE, GIVEN THE NUMBER OF PULSES n, THE PROBABILITY OF 
FALSE ALARM 
% PFA AND THE AVERAGE SIGNAL TO NOISE RATIO S/N. 
%* EXAMPLE: SWERLS5(PFA,N,SBAR) ; 
Eke ekee eee eee eee eee eee Kee eee ere reek ee eee Ah 
Emmet Lon pd=Sswerls(pfa,n,sbar) 
Yb=thresh(pfa,n) ; 
Z=n*sbar; 
YM=exp (-yb) ; 
OK = Q; 
r= } 
YMS=YM; 
fer 1=1:k-1 
YM=YM*yb/i; 
YMS = YMS + YM; 
end 
SUM=YMS ; 
XB=exp(-Z) ; 
XBS=XB; 
while OK == 0 
YM=YM*yb/k; 
YMS = YMS + YM; 
SUM= SUM+ YM*(1-XBS) ; 
el=1-YMS; 
XB =XB*2z/(k-n+1) ; 
XBS = XBS + XB; 
e2=1-XBS;- 
er=el*e2; 
if er <= le-6 
OK=1; 
end 
k = k + 1; 
end 
pa = SUM; 
return; 


% FILE NAME: THRESH.M 

Brewer ekaekekkekeeehkekkekkke keke ete eee Keer k eee eee eK eke ee ee ee 

%** RETURNS THE THRESHOLD yb FOR A GIVEN PROBABILITY OF FALSE 
ALARM (pfa) 

%* AND A GIVEN NUMBER OF PULSES 

%* EXAMPLE: THRESH(Pf£a,N). 

Eee ekeeew kee eee eee kere keke keke eee eee ere eee eee ee ee a we we 

function y=thresh(pfa,n) 

format long 


meen =< 1 | n==fix(n), 
error('Wrong input number of pulses'); 
break 

end 

ee = 1 
Y==leg(pta); 


a1 


recur 
end 
l = =logl0(pEa); 
N2=Sore( nj lease eae, 
Ven-n24+2 35 22 el oe 
comp=le-6; 


Yatio= |: 
while comp <= ratio 
p=prob(n,y) ; 
ym=1; 
for kai ne te 
ym=ym/k*y ; 
end; 
ym=ym* exp (~y) ; 
dely=(p/ym) *log(p/pfa) ; *%* correction magitude 
y=y+dely; 
ratio=abs(dely/y) ; 
end 
return 


%* FILE NAME: PROB.M 


EeNKK Kee kee eek eee eee ee 


5 THIS PROGRAM RETURNS THE FALSE ALARM PROBABILITY P(n, yb) 
as in _ 
%$*  SUMMATORY (yb*k/fact(k)) for k=0 to n-1, AND IS A COMMON 


FUNCTION sj 
%** FOR THE PROBABILITY OF DETECTION COMPUTATIONS 
” 


5 Ee SaiP EE PROB (iN, xe) 
* 


Beet kee eek e ree eeerk ere e ree ek eee e ee keke ee eree  e e 


EUNGCELOM p=OroD (x,y) 
PE oO ese — (oe), 
error('Number of pulses should be integer and greater than 
zero'), break, end 
Htf==1), 79769313 4862069 =lossos- 
E=709.782712893384040999999; 
p=0; 
2] <= 
return; 
else: £ x== 
p = exp(-y); 
return; 
else 
%* -[1£ x>1] 
t=l1; 
for k=): 1-x-1; 
ESt/ key: 
P=C+p; 
end 
p=(p+1) *exp(-y) ; 
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=~: Te eer 


end 
return; 


*File name: THRESHM.M 
function yb=threshm(pfa,n) 
tol=le-6; 
E=709.782712893384040999999; 
1 = -logl0O(pfa) ; 
me—sqrt (mn) ;i2Z=sqrt(l); 
yen-n2+2.3*12* (12+n2-1) ; 


Plog=logl10(pfa); 
xpfa= 1/pfa; 
ok=0; 

while ok==0 


me y <= & 


m= 0 
yix=0; 
ylog=log(y) ; 
tsum=exp(-y) ; 
EOram —=.0-13n-1 
m=m+1 ; 
sunk=log(m) ; 
ylx= ylx+ylog -sunk; 
tsum=tsum+exp (ylx-y) ; 
end 
yl=y+log(xpfa*tsum) *tsum/exp(ylx-y) ; 


if abs(y-yl) < tol 
yb=yl; 
ok=1; 
return 
else 
y=yl; 
ok=0; 
end 


else 
m=0- 
yilx=0; 
ylog=log(y) ; 
okk=0; 
while okk==0 


m=m+1; 


sunk=log(m) ; 
ylx=ylx+ylog-sunk; 
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1£ y-vylx <= 5 


tsum=exp(ylx-y) ; 
kk — 1 
else 
okk=0; 
end 
end 
ee i oo 
yl=yt+log(xpfa*tsum)* tsum/exp(-y+ylx) ; 
if abs(yl-y) < tol 
yb=y1 ; 
Oka 
break; 
return 
else 


else 

for m=m:l:n-l 
sunk=log(m) ; 
ylx=ylx+ylog-sunk; 
tsum=tsum+exp (-y+ylx) ; 

end 

yl=y+log(xpfa*tsum) *tsum/exp ( “y+yl1x) ; 

if abs(yl-y) <= tol 

yb=yl; 
ok=1; 
break; 
return 


end 
end 
end 


*FILE NAME: MARCUM.M 
EUNGCELOM LD, vmo] —ocoumalian..) 
format long 

err=le-6; 
E=709.782712893384040999999; 


%*eeeee Start the program **tttrert renee eee ee eee rene ee 
lamdal=(n/(2*y) )*2+(n*x) /y; 


lamda=1-(n/(2*y) )-sqrt(lamdal) ; 
auxl=(-lamda*y)+(n*x*lamda) /(1-lamda) ; 
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auxZ=exp (auxl) /(1-lamda) “n; 
1f aux2 <= err 
aux3=n* (x+1)>; 
1f y >= aux3 
p=0;return 


else 
p=l;return 
end 
else 
lb=n-1; 
k=0; 
x=n*x; 
Me x <= £ 
xk=exp (-x) ; 
xks=xk; 
m=O - 
1€ y<= E 
=exp(-y); 
yms=ym,; 
SU a ara ea a ah as eevee ee eae a ae we 


for m=l:n-1l 
ym=ym*y/m; 
yms=yms+ym; 
a err>(l-yms) 
Pp=yms ; 
ymo=ym 
x=x/n; 


return 
end 
end 
m=n-1; 
SUM=yMs ; 
ymo=ym 
ok=0; 
while ok==0 
k=k+1; 
m=m+1 ; 
ym=ym*y/m; 
yMs =yMS+ym ; 
sum=sum+ym* (1-xks) ; 
xk=xk*x/k; 
xks=xks+xk; 
ans=(l-yms) * (1-xks) ; 
1f ans<= err 
ok=1; 
end 
end 


a) 


p=sum; ymo=ym; 
else 
ymly=0; 
Viog=-loagi)] 
okk=0: 
while okk== 
m=m+ 1 ; 
Smt = | Ogam).- 
ymly=ymly+ylog-xmf; 
1£ (y-ymly) <= E 
ym=exp (-yt+ymly) ; 
OkK=a, 
else 
okk=0; 
end 
end 
1£ om <= 915 
yms=ym; 


1£ err>(l-yms) 
Pp=yms ; 
ymo=ym 
x=x/n; 
ook=1; 
return 

end 


end 
SUM=yms ; 
ymo=ym 
Ok=0- 
while ok==0; 
k=k+1; 
m=m+1; 
ym=ym*y/m; 
yMms=yms+ym; 
sum=sum+ym* (1-xrs) ; 
xk=xk*x/k; 
XrS=xXrs+xk; 
ans=(l-yms) *(1-xrs) ; 
1f ans > err 
ok=0; 
else 
p=sum; 
x=x/n; 
break; 
return; 
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end 
p=sum; ymo=ym; 


else $------------ fOr m>i1b; 
ymo=0; return 
yms=ym; 
sum=0; 
kd=m-n; 
Die ke 
p=0; 
x=x/n; 
break; 
return; 
else 
okkk=0; 
while okkk==0; 
k=k+1; 
kk=xk*x/k; 


xks=xks+xk; 
while (k-kd) < 0 
okkk=0; 
end 
end 
ooook=0; 
while ooook==0; 
sum=sum+ym* (1-xrs) ; 
xk=xk*x/k; 
xrs=xrs+xk; 
ans=(1l-yms) *(1-xrs) ; 
while ans < err 
p=sum; 
x=x/n; 
ooooke=l1; 
break; 
return 
end 
k=k+1; 
m=m+1 ; 
ym=ym*y/m; 
YMS=yms+ym; 
oo00ok=0; 
end 


end 
end 
end 


else 
x<kix=0- 


od 


xlog=log(x); 
o000k=0; 
while oook== 
k=k+1; 
xklx==xklx+xlog-xkf€; 
while (x-xklx) > E 
oook=0; 
end 
end 
xk=exp (-x+xk1lx) ; 
1b=ib+k; 


1£f y<= E 
ym=exp(-y); 
yms=ym,; 


yms=yms+ym; ° 
1£ err>(l-yms) 
P=yMs ; 
ymo=ym; 
x=x/n; 
return 
end 
end 
SUM=ymMs ; 
ymo=ym; 
ok=0; 
while ok==0; 
k=k+1; 
m=m+1; 
ym=ym*y/m; 
yms=yms+ym; 
sum=sum+ym* (1l-xrs) ; 
MK=KK "xk: 
XOS=Xrs+xk; 
ans=(l-yms) *(1-xrs); 
1f ans<= err 
ok=1; 
end 
end 
p=sum; 


else 
ymly=0 ; 
ylog=log(y); 
okk=0; 
while okk== 
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m=m+ 1 ; 
ele = Oat). 
ymly=ymly+ylog-xmf; 
1f (y-ymly) <= E 
ym=exp (-y+ymly) ; 


okk=1; 
else 
okk=0; 
end 
end 
ie fete oe 
yms=ym; 
ook=0; 


while ook==0 
while m~=1b 
m=m+1 ; 
ym=ym*y/m; 
YMS=yMms +yM,; 
1f err>(l-yms) 
D=ymMs ; 
ymo=ym; 
x=x/n; 
ook=1; 
break; 
return 
else 
ook=0; 
end 
end 
end 
SUM=yMs ; 
ymo=ym; 
ok=0; 
while ok==0; 
k=k+1; 
=m+1; 
ym=ym*y/m; 
YMS=yms+ym; 
sum=sum+ym* (1-xrs) ; 
xk=xkK*x/k: 
XOS=Xrs+xk; 
ans=(l-yms) *(1-xrs) ; 
1f ans > err 
ok=0; 
else 
p=sum; 
x=x/n; 
break; 
return; 


a9 


end 
end 
else %------- +--+ - EOL. mao 
ymo=0 ; 
yms=ym; 
sum=0 ; 
kd=m-n; 
if kd>x 
p=0; 
x=x/n; 
break; 
return; 
else 
okkk=0; 
while okkk==0; 
k=k+1; 
Kk=xk tx/ kK 
xks=xks+xk; 
while (k-kd) < 0 
okkk=0; 
end 
end 
COOOK=E; 
while ooook==0; 
sum=sum+ym* (1-xrs) ; 
xk=exk*x/k; 
XYS=xXrst+xk; 
ans=(l-yms) *(1l-xrs) ; 
while ans < err 
p=sum; 
x=x/n; 
OOCeK= 1; 
break; 
return 
end 
k=k+1; 
m=m+1 ; 
ym=ym*y/m; 
yms=yms+ym; 
oo0ook=0; 
end 
end 
end 
end 
end 
end 
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